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G^ ■ Abstract 

^^ I We propose a modification of the Hybrid-Monte-Carlo algorithm that 

^— V ' allows for a larger step-size of the integration scheme at constant accep- 

tance rate. The key ingredient is that the pseudo-fermion action is split 
^^ I into two parts. We test our proposal at the example of the two-dimensional 

"t^ ' lattice Schwinger model with two degenerate flavours of Wilson- fermions. 
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1 Introduction 

The Hybrid-Monte-Carlo algorithm |l|] has become the standard algorithm to 
simulate lattice QCD with dynamical fermions. E.g. it has been used in the 
recent large scale study of the spectrum of light hadrons. 

At present, most simulations are performed for two flavours of mass-degenerate 
sea-quarks, where the mass is by far larger than the masses of the up and down 
quark. An elaborate extrapolation of the data towards the chiral limit is needed. 

Therefore it is desirable to further reduce the mass of the sea-quarks. However, 
reducing the mass of the sea-quarks, i.e. approaching n^. in the case of Wilson 
fermions, the cost of the simulation increases for at least three reasons (see e.g. 
ref. i): 

• The condition number of the Dirac-matrix increases. Hence the number of 
iterations that are needed for the inversion of the Dirac-matrix increases. 

• The autocorrelation times in units of trajectories increase as Kc is ap- 
proached. 

• The step-size of the integration scheme has to be decreased as k^ is ap- 
proached to maintain a constant acceptance rate. This means that for one 
trajectory more inversions of the Dirac matrix have to be performed. This 
effect can be nicely seen in table II of ref. . 

While much work has been devoted to deal with the first two problems, little 
attention has been paid to the third. Recently it has been noticed that the step- 
size at a given acceptance rate depends on the precise form of the pseudo-fermion 
action @, ^. It was shown that preconditioning of the fermion matrix leads to 
larger step-sizes in the Hybrid-Monte-Carlo algorithm. 

In this paper we propose a rather simple modification of the pseudo-fermion 
action that allows to increase the step-size at fixed acceptance rate. In particular, 
we split the pseudo-fermion action into two parts, separating (partially) the small 
and the large eigenvalues of the Dirac-matrix. 

The paper is organised as follows. In section 2 we introduce the modi- 
fied pseudo-fermion action. We test our proposal at the example of the two- 
dimensional two-flavour Schwinger model. In particular, we show that our idea 
can be applied in addition to even-odd preconditioning that is discussed in section 
3. In section 4 we present our numerical results. Finally we give our conclusions 
and an outlook. 



2 The modified pseudo-fermion action 

The following discussion is rather general and applies, e.g., to four- dimensional 
lattice QCD with Wilson fermions as well as the two-dimensional Schwinger 
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model on the lattice. We start our discussion with the partition function for 
two degenerate flavours of Wilson fermions: 

Z= fD[U]expi-SGiU)) detM(f/)^M(f/) , (1) 

where U is the gauge field, Sg{U) the gauge action and 

M{U) = 1 - kH{U) (2) 

the fermion matrix, where k is the hopping parameter. 

In the Hybrid- Monte-Carlo algorithm, the determinant detM^M is repre- 
sented by the integral over an auxiliary field (pseudo-fermions) (p: 

detM^M oc /"D[0]D[0"f] exp(-|M-Vr) • (3) 

Hence the action, as a function of the gauge-field and the pseudo-fermion fields, 
is given by 

^(t/,0) = SciU) + SpiU.ct)) (4) 

with 

^^(t/,0) = iM(f/)-vr • (5) 

Our new idea is to split the fermion matrix into two factors. Each factor is 
represented by an integral over an auxiliary field: 



(6) 



detM^M oc /D[^]D[7/'t] D[0]D[(/)^] exp(-|M^Vn exp(-|MM-Vn 

where the auxiliary fermion matrix M is chosen as 

M = 1-RH (7) 

with K < K. I.e. the pseudo-fermion action is now given by the sum of the two 
terms 

5^i(f/,^) = |M(f/)-Vl' (8) 

and 

SF2iU,^) = |M([/)M([/)-V|' . (9) 

The condition number of M as well as i? = MM~^ is reduced compared with the 
original matrix M. This is the reason, why we expect, similar to preconditioning 
[^ O, that the step-size in the Hybrid-Monte- Carlo can be increased. In the case 
of M the condition number is reduced for typical gauge-configurations since we 
have chosen k < n. Next we consider R. Let us write 

M = aM + b (10) 



with 

a = — , 6=1— a . (11) 

K 

Hence 

R~^ =a + bM-^ . (12) 

For K close to Kc, lAmmT^ > 1 and \Xmax\ = 0(1), where Xmin and Xmax are 
the minimal and maximal eigenvalues of M. Hence the condition number of R is 
essentially reduced by a factor of b compared with the condition number of M. 

In the Hybrid-Monte-Carlo, the variation of the action with respect to the 
gauge field has to be computed. In the case of the standard pseudo-fermion 
action (||) one obtains 

SSf = -0^ [M^-'^Ar'^ 6M M~^ + M^~^ 5M^ M^-^M'^] (/. 

= -\Y^ 5M X + X^ 5M^ Y] , (13) 

where X = M-^(j) and Y = M^~^X. 

In the case of the split pseudo-fermion action 6Sfi can be computed exactly 
as in eq. (0). Also 6Sf2 can easily be computed. 

SSf2 = -<P^ [R^-^R-^ 6R R-^ + R^-^ 6R^ R^-^R-^] </) 

= -0"^ [{a + bM^-^)M-^ 6M M"^ + M^"^ 5M^ M^~\a + bM'^)] 
= -[yUM X + X^ 5M^ Y] , (14) 

where we have used R~^ 6R R~^ = b M~^ SM M~^ which follows from eq. (p!^. 
The auxiliary vectors X and Y are modified compared with eq. (plsl) : X = M~^0 
andy = Mt-i(a0 + 6X). 

I.e. the computational cost to compute the variation of the second part of 
the modified pseudo-fermion action (P) is much the same as for the standard 
pseudo-fermion action (|^). In both cases we have to apply M~^ and M^~^ to a 
vector. 

Computing the variation of the first part of the modified pseudo-fermion ac- 
tion (H) means an overhead compared with the standard case. However, since 
k < K fewer iterations are needed to compute M~^ip than M~^(f). 

We also should note that the modification can be very easily implemented, 
given a standard Hybrid-Monte- Carlo program is at hand. One just has to add 
a second pseudo-fermion field for R and implement the modified definition of the 
auxiliary vectors X and Y. 

In the following simulations we like to test whether the modified pseudo- 
fermion action indeed allows for a larger step-size at constant acceptance rate 
and whether this larger step-size outweighs the computational overhead discussed 
above. 



3 The Schwinger model 

We have tested our proposal to split the ferniion matrix in the two dimensional 
Schwinger model with two degenerate flavours of Wilson fermions. It has been 
used frequently as a toy-model to study properties of Monte Carlo algorithms. 
The gauge action of the Schwinger model is given by 

Sg = -/?^Re Upiaq,x , (15) 

X 

where 

Uplag,x = Ux,0 Ux+(lfl),l f^:^+(0,l),0 Ux,l 5 (16) 

where the link variables Ux^^ are elements of U{1) and U* is the complex con- 
jugate of Ux^fM- X labels the points of the two-dimensional square lattice and /i 
gives the direction. The lattice constant is set to a = 1. The fermion matrix can 
be written as 

M = 1-kH , (17) 

where the hopping part of the fermion matrix is given by 

H = Y,{S.-^,ya + l,)Ux-f.,, + 5x+^,ya-l,)U:,,) , (18) 

where ft is unit vector in /i direction. The 7-matrices in two dimensions are given 
by the Pauli-matrices a. In our simulations we applied even-odd preconditioning. 
See e.g. ref. [Q. The sites of the lattice are decomposed in even and odd sites. 
Then the fermion matrix can be written in the form 

where H^o connects odd with even sites and Hoe vice versa. For the fermion 
determinant the identity 

detM = det(lee - n^HeoHoe) (20) 

holds. Hence the original problem is reduced by half in the dimension (of the 
fermion matrix). The pseudo-fermion field which is used for the stochastic 
estimate of the fermion determinant lives only on even sites. The even-odd pre- 
conditioned fermion matrix is given by 

Mee = lee-t^^ Heo Hoe • (21) 

In the following we shall apply our proposed splitting (^) to the even-odd pre- 
conditioned fermion matrix Mee- We use 

Mee = lee — i^ Heo Hoe ■ (22) 

I.e. in eq. (|TT]) k, and k have to be replaced by k^ and k"^. 



4 Numerical results 

In our simulations with the Hybrid-Monte-Carlo algorithm we have used the 
standard leap-frog integration scheme. In order to separate the effects of the 
gauge action and the pseudo-fermion action we applied the two step-size method 
of Sexton and Weingarten p. In particular, we have chosen n = 4 as ratio of 
the step sizes. The length of the trajectory was chosen randomly between 0.5 
and 1.5. We computed M'^cf) with the BiCGstab algorithm. In the rare cases, 
where BiCGstab did not converge, we used the conjugate gradient algorithm. We 
started the iteration with (p. As stopping criterion, we required that the residue 
is less that 10"^°. 

We simulated at parameters (3 and k that were recently used by Peardon [Q 
to study the effects of preconditioning on the performance of the algorithm. 

We performed a first set of simulations for L = 32, /3 = 4.0 and k, = 0.26. 
In table |l] we have summarised results of simulations using various values k. 
In all these simulations we generated 20000 trajectories. 1000 trajectories were 
discarded for thermalisation. In all our simulations we fixed the acceptance rate 
to about 0.8. For this purpose we performed a few preliminary simulations for 
each K with smaller statistics. 

For k = 0.0 we obtained an acceptance rate of about 0.8 with a step-size of 
dr = 0.035. Here dr is the step-size used for the pseudo-fermion action. The 
step-size for the gauge-action is dr/4. Increasing k the step-size can be increased 
up to about dr = 0.07 for k = 0.22. Increasing k to 0.255 leads to a smaller step- 
size again. This increase is not surprising, since in the limit k ^ k we recover 
the standard algorithm again. 

The number of iterations m and m needed to compute M~^%p and M~J^(j), 
respectively, refer to the inversions performed for the accept/reject step at the 
end of the trajectory. I.e. they are not performed for equilibrium configurations. 
Therefore it is not surprising that m depends slightly on k. 

For k = 0.22 the numerical overhead to compute the modified fermion action 
compared with the standard one is rather moderate. The number of iterations 
needed to compute M'^ip is only a small fraction of the iterations needed to 
compute Mg~V. 

We have measured square Wilson-loops up to size of 5 x 5 and the topological 
charge, using the geometric definition. The estimates of the expectation values 
were consistent among the simulations that we have performed. Also, the results 
for the Wilson-loops of size 1x1 and 4x4 are consistent with those quoted by 
Peardon Q] in his table 2. It turned out that the auto-correlation times (in units of 
trajectories) of the various observables are, within statistical errors, independent 
of k. Therefore, the runs can be compared solely on the ground of the numerical 
cost per trajectory. 

The numerical costs are dominated by applications of M^e and Mjg. Therefore, 
as a measure of the computational cost we took the number of applications of M^e 



Table 1: Results for L = 32, f3 = 4.0 and k, = 0.26. In the first column we 
give the hopping parameter R of the auxiliary fermion matrix M^e ■ The second 
column contains the step-size dr. In the third and fourth column, the number 
of iterations needed to solve M~J-(f) and M~^ip is given, respectively. In the fifth 
column we give the acceptance rate. Finally, in the sixed column the cost, which 
is defined in the text, is presented. Details are discussed in the text. 



K 


dr 


m 


m 


ace 


eost 


0.0 


0.035 




203.42(21) 


0.808(4) 


24788(53) 


0.18 


0.056 


14.576(5) 


204.85(16) 


0.793(4) 


17107(35) 


0.20 


0.063 


19.771(5) 


205.42(13) 


0.803(3) 


15724(31) 


0.21 


0.067 


23.724(6) 


205.59(13) 


0.784(4) 


15120(29) 


0.22 


0.07 


29.440(8) 


206.08(13) 


0.784(4) 


14898(29) 


0.23 


0.072 


38.394(11) 


206.26(13) 


0.784(3) 


15058(28) 


0.24 


0.07 


54.596(21) 


206.61(12) 


0.801(3) 


16515(33) 


0.25 


0.07 


91.91(6) 


206.84(16) 


0.794(3) 


18958(37) 


0.255 


0.065 


134.29(10) 


207.09(12) 


0.774(3) 


23181(47) 



Table 2: Results for L = 64, (3 = 4.0 for the two values of k = 0.257 and 
K = 0.2605. For k, = 0.257 we have performed 10000 trajectories for each value 
of K. For K = 0.2605 we have performed 6600 trajectories for k = 0.0 and 8400 
trajectories for k = 0.23. The notation is the same as in table 1. 



K 


R 


dr 


m 


m 


ace 


cost 


0.257 
0.257 


0.0 
0.22 


0.04 
0.054 


30.210(10) 


195.8(2) 
202.4(2) 


0.784(5) 
0.783(5) 


21669(71) 
18681(55) 


0.2605 
0.2605 


0.0 
0.23 


0.022 
0.051 


39.69(2) 


384.4(8) 
398.2(7) 


0.792(7) 

0.788(5) 


74550(330) 
37038(130) 



and Mjg per trajectory. The minimum of costs is rather shallow and is located 
around k = 0.22. Compared with the standard simulation {k = 0) we see a 
reduction of the cost by a factor of 24788/14898 ^ 1.66. 

Finally we performed simulations for L = 64, j3 = 4.0 at k = 0.257 and 
K = 0.2605. Following Peardon ^, the pseudo-meson masses at these values of 
K are am,p = 0.210(3) and amp = 0.124(5), respectively. 

Here we performed only simulations for the standard case k = and for one 
non-trivial value of k. We have chosen k = 0.22 for k, = 0.257 and k = 0.23 for 
K, = 0.2605. The simulations for the L = 32 lattice have shown that the minimum 
in the computational cost is rather shallow in k. Therefore, we expect that the 
computational cost at our values of k is rather close to the optimum. 

First of all we notice that for k = the step-size dr has to be decreased when 



K is increased. For k 7^ the decrease of the step-size with increasing k is much 
smaller. I.e. the performance gain that we achieve with the modified pseudo- 
fermion action increases as Kc is approached. For k = 0.2605 with k = 0.23 we 
get an improvement of a factor 74550/37038 ~ 2. 

5 Conclusions and Outlook 

We have proposed to use a modified pseudo-fermion action in the Hybrid Monte 
Carlo simulation. This modification is based on a factorisation of the fermion 
matrix. The modification can be easily implemented, given a code for the stan- 
dard Hybrid-Monte- Carlo algorithm is at hand. We have tested our proposal at 
the example of the two-dimensional two-flavour Schwinger model with Wilson 
fermions. The numerical tests have shown that the modiflcation of the fermion 
action indeed allows for a larger step-size of the leap-frog integration scheme at a 
fixed acceptance rate. Moreover, the increased step size outweighs the numerical 
overhead due to the modified action. For our largest value of k, we found a net 
gain in performance of a factor of two. 

There remain a number of open problems. First of all we would like to test 
our proposal for lattice QCD. We would like to study more systematically the 
dependence of the improvement on the parameters (3, k and the lattice size. Also 
a better theoretical understanding of how the performance of the Hybrid-Monte- 
Carlo depends on the form of the pseudo-fermion action would be helpful. It 
would be interesting to test how higher-order integration schemes [|] perform in 
combination with the modified pseudo-fermion action. Also, one might incorpo- 
rate our new idea in the Polynomial- Hybrid- Monte-Carlo algorithm [^ §. 
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